Why twisting angles are diverse in graphene Moire patterns? 
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The interlayer energy of the twisting bilayer graphene is investigated by the molecular mechanics 
method using both the registry-dependent potential and the Lennard- Jones potential. Both po- 
tentials show that the interlayer energy is independent of the twisting angle 8, except in the two 
boundary regions 6 ~ or 60°, where the interlayer energy is proportional to the square of the 
twisting arc length. The calculation results are successfully interpreted by a single atom model. An 
important information from our findings is that, from the energy point of view, there is no prefer- 
ence for the twisting angle in the experimental bilayer graphene samples, which actually explains 
the diverse twisting angles in the experiment. 

PACS numbers: 61.48.Gh, 61.72.Nn 



I. INTRODUCTION 

The electronic band structure in single-layer graphene 
is well-known for its linear Dirac cones around K points 
in the Brillouin zone. This Dirac electron can be well de- 
scribed by a single 7r-orbital tight-binding model- Differ- 
ent from the single-layer graphene, a parabolic electronic 
dispersion was predicted theoretically for the Bernal- 
stacking bilayer graphene (BLG)»i— However, several ex- 
perimental groups have observed the Dirac-like electron 
in BLG through direct or indirect measurements, which 
is attributed to the twisting defect in the BLG4~— The 
twisting pattern can be observed by the scanning tunnel- 
ing microscop y 12 ! 14 ' 15 , high-resolution transmission elec- 
tron microscop y 16 ' 17 , and atomic force microscopy J£ Be- 
sides electronic properties, it has been widely shown that 
single-layer graphene and BLG have peculiar mechanical 
and lattice properties J^r— There is increasing interest 
in studying possible effects from the twisting pattern on 
the lattice properties of the BLG i 18 i 25 ~ — The most re- 
cent experiment demonstrates that the Raman spectrum 
strongly depends on the twisting angle of the Moire pat- 
tern in BLG^ 6 - 

In above studies of various properties of the twisting 
BLG, a fundamental issue is to prepare the twisting an- 
gle of the BLG sample. The observed twisting angles are 
diverse in existing experiments. For example, in Ref. @, 
Luican et.al obtained a twisting angle around 21.8° or 
3.5°. In Ref. [T3, Brown et.al have observed a broad 
distribution of the twisting angle around 29°, 24°, 17°, 
12°, and 5°. The twisting angle is measured to be 4° 
in Ref. [lfjl Now, a fundamental question arises: Is there 
any preference for the twisting angle in the BLG sample? 
Should the experiment always observes a commensurable 
angle in the twisting BLG? The present work studies the 
interlayer interaction in the twisting BLG to shed some 
light on these questions from the energy point of view. 

It has been a great challenge for theoretician to calcu- 
late the interlayer interaction in layered structures like 
BLG. To calculate the atomic energy, one may apply 



cither density functional theory (DFT) or an empirical 
potential. In the empirical potential, the van der Waals 
(vdW) interaction in such layered structure is usually de- 
scribed by a Lennard- Jones potential. However, on the 
one hand, the interlayer interaction is long-range, so the 
standard DFT approach can not describe it, because the 
DFT is based on a local density approximation or a gen- 
eralized gradient density approximation^ On the other 
hand, it has been pointed out that the vdW potential 
itself is not sufficient to describe the interlayer energy, 
especially for the interlayer shearing movement, because 
the shearing is dominant by the 7r-overlap between differ- 
ent layer s 30 ' 31 and the registry matching plays an impor- 
tant role^ 2 - The shearing property calculated from a pure 
vdW potential is one order smaller than the experimental 
value in such layered structure! 33 ' 34 

There are mainly two solutions for this issue. The 
first method is to develop a vdW-corrected DFT (DFT- 
D) approach, where the long-range interaction is de- 
scribed directly by the vdW potential 3 - 5 - - — or is in- 
cluded through a density-density interaction in the DFT 
scheme 4 ^. The second method is to include the tt- 
overlap through some empirical potential terms with 
empirical parameters fitted to experiment or DFT-D 
results i^r—iiir— As pointed out by Girifalco and Hodak 
in 2002, these two methods are actually related to each 
other and should give the same results for interlayer in- 
teraction if they are applied properly^ In present study, 
the registry-dependent empirical interlayer potential will 
be applied in the calculation of the interlayer interaction, 
as long as we aim at simulations for large systems. 

In this paper, we calculate the interlayer energy in the 
twisting BLG. The interlayer interaction is described by 
cither the registry-dependent potential or the Lennard- 
Jones potential. We find that the interlayer energy does 
not depend on the twisting angle, except in the boundary 
regions 9 m 0, 60°, where the interlayer energy is pro- 
portional to the square of the twisting-related arc length 
S = R9. We explain these results by averaging the energy 
distribution in a single atom (SA) model. The twisting 
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FIG. 1: (Color online) The interlayer potential energy calcu- 
lated from the registry-dependent potential with different cut- 
off. The y-axis (E a ) is the total interlayer energy per atom. 
In this calculation, the radius of the top graphene layer is 
R = 100 A. The twisting angle 9 = 0.0, i.e AB-stacking BLG. 
For a cut-off r c = 25 A, the variation in the energy is on 
the order of 10 -3 meV. Inset displays the cross section of a 
twisting BLG, where R is the radius and r c is the boundary 
region. 

anglc-dcpcndcnce of the interlayer energy actually pro- 
vides an explanation for the diverse twisting angles ob- 
served in the experiment. 

The present paper is organized as follows. In Sec. II, 
after a brief description of the structure, we present the 
detailed results based on the registry-dependent poten- 
tial. Sec. Ill is devoted to the simulation results for the 
Lennard- Jones potential. The paper ends with a brief 
summary in Sec. IV. 



II. RESULTS AND DISCUSSION FOR 
REGISTRY-DEPENDENT POTENTIAL 

We start with the AB-stacking BLG. The xy axes sit in 
the bottom layer. In the AB-stacking, there is a carbon 
atom from the top layer sitting on the head of a carbon 
atom from the bottom layer. The z axis comes across this 
pair of atoms, and the origin of the coordinate system is 
set on the carbon atom in the bottom layer. The bot- 
tom layer is infinite large. This 'infinite' is numerically 
realized by ensuring that the radius of the bottom layer 
is always larger than (i? + r c ), where R is the radius of 
the top graphene layer and r c is the cut-off in the inter- 
layer potential (see inset in Fig.QJ. The top layer is then 
twisted about the z axis for an arbitrary twisting angle 9. 
In particular, the AB-stacking BLG has a twisting angle 
9 = 0, and the AA-stacking BLG has 9 = 60°. 

Previous research has shown the limitation of the 
Lennard-Jones potential in the description of the inter- 
layer interaction in BLG^ Particularly, it only gives less 
than 10% of the twisting or shearing energy in BLG. 



FIG. 2: (Color online) The saturation of the energy per atom 
with increasing radius for different twisting angles 9. (a) 
9 = 0, i.e AB-stacking BLG. Solid and open dots correspond 
to cut-off r c = 25 and 50 A. (b) 9 = 60°, i.e AA-stacking 
BLG. The energy per atom in the AA-stacking BLG is radius 
independent, (c) 9 = 21.8°, i.e the first commensurable angle, 
(d) 9 = 13.2°, i.e the second commensurable angle. 

To improve the situation, the Lennard- Jones potential is 
extended to be registry-dependent £2r— In our calcula- 
tion, we have employed the latest version of the registry- 
dependent potential developed by Lebedeva et.alin 2011, 
which has been succeeded in predicting the energy bar- 
rier for a relative translation of the two graphene layers 
in BLGi 34 ' 46 In this model, the interaction between two 
atoms from adjacent graphene layers is described by fol- 
lowing formula: 

V(r) = A^~y + Be- a( - r - z ^ 

+ C (1 + D lP 2 + D.y) e- x ^ 2 e- x < z2 - z o),(l) 

where r is the distance between two atoms, and p 2 = r 2 — 
z 2 . The parameters are as follows: A = —10.510 meV, 
z Q = 3.34 A, B = 11.652 meV, a = 4.16 A" 1 , C = 
35.883 meV, D Y = -0.86232 A" 2 , D 2 = 0.10049 A" 4 , 
Ai = 0.48703 A~ 2 , A 2 = 0.46445 A" 2 . The interaction 
cut off is r cut = 25 A. 

For different twisting angle, the structure is relaxed to 
the energy minimum state with boundary atoms fixed. 
In particular, the interlayer space is also relaxed. We 
then calculate the interlayer energy per atom in this op- 
timized structure, E a = £ , to t/A ? top/2, where -Etot is the 
total interlayer energy between the two graphene layers 
and A^top is the total carbon atoms in the top layer. A 
factor of 2 is to divide the pair energy into each atom. 
We note that E a here is different from the binding en- 
ergy defined in Refs. [34U46I by a factor of 2. We first 
examine effects from the cut off in the potential. Fig. [1] 
shows the energy per atom calculated with different cut 
off in AB-stacking BLG. In this calculation, the radius 
of the top layer is 100 A. The variation is on the order 
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FIG. 3: (Color online) The energy per atom versus twisting 
angle 9. The whole curve is symmetric about 9 = 60° due to 
the mirror symmetry between the two nonequivalent atoms 
in graphene, so 9 £ [60, 120]° is not shown. Inset shows the 
interlayer space versing twisting angle 9. 

of 10~ 3 meV for cut-off 25 A. In the following calcula- 
tion for the registry-dependent potential, we use a cut-off 
25 A, so that our calculation precision is 10~ 3 meV. 

Fig. [5] shows the convergence of E a with increasing ra- 
dius (R) in BLG with various twisting angles. The twist- 
ing angles are 0.0, 60°, 21.8°, and 13.2° in panels from (a) 
to (d). These angles correspond to the commensurable 
angles from Santos's equation^ 8 = arccos((3i 2 + 3i + 
0.5)/(3z 2 + 3i + 1)) with integers i =0, 1, 2, and +oo. 
For all panels, the energy saturates at radius around 
100 A, i.e the oscillation in the energy for a system 
with r > 100 A is smaller than our calculation preci- 
sion (10 -3 meV). In panel (a), i.e AB-stacking BLG, E a 
converges to a saturate value of -23.345 meV. In this cal- 
culation, two different cut-off have been used. The data 
corresponding to r c = 25 A are represented by solid dots. 
Results for r c = 50 A are displayed by open symbols. 
The saturate value is the same in both case. It further 
confirms that r c = 25 A is suitable. In panel (b), i.e 
AA-stacking BLG, E a =-17.215 meV is radius indepen- 
dent. The energy difference between these two stacking 
style is 6.130 meV. In panels (c) and (d), E a converges to 
the same saturate values of -20.654 meV. It is interesting 
that the energy in a BLG with different commensurable 
angles are the same. It indicates that the commensu- 
rable twisting angle may not correspond to the energy 
minimum state. Hence, from the energy point of view, 
the twisting angle existing in the experiment is not nec- 
essarily the commensurable angle. 

Fig. [3] further confirms that E a does not depend on 
the twisting angle in a wide angle range 8 £ [10,50]°, 
except in the two boundary regions (8 rj 0, 60°). For 
small radius, there is some obvious oscillation in the 
range of 8 € [10,50]°. With increasing radius, this os- 
cillation amplitude decays and becomes indistinguishable 
after R > 50 A. It indicates that the oscillation is actually 



-21 




0.2 0.4 0.6 

Twisting angle (degree) 



FIG. 4: (Color online) The energy per atom versus twisting 
angle around 9 — 0. For curves from bottom to top, the radius 
increases from 10 to 100 A. All curves can be well fitted to 
E a = E° + a9 2 , where E° is the value at 6 = 0. a is the 
fitting parameter. Inset shows the radius-dependence of a. a 
is fitted to a = 9.9 x 10~ 4 R 2 . 



due to the size effect, and is not related to the twisting 
angle. Similar size effect also exists in a recent work by 
Shibuta and Elliott 42 The inset in Fig. [3] shows the opti- 
mized interlayer space of the BLG with different twisting 
angle, which looks quite similar as the energy curve. Our 
calculation predicts that the interlayer space is not sen- 
sitive to the twisting angle in a large angle range. 

In Fig. [31 although the curve is a platform in a wide 
angle range, we can see that it changes sharply around 
the two boundaries 9 rs and 60°. These two regions 
are zoomed in in Fig. 0] and Fig. [5] Fig. [4] shows that E a 
can be fitted to E a = E® + a9 2 , where E® is the value at 
8 = 0. a is the fitting parameter and is radius-dependent. 
Inset shows the radius-dependence of a, which is fitted 
to a = 9.9 x 10~ 4 i? 2 . Hence, around 8 = 0, we have 
E a = E Q a + 9.9 x 10- 4 i? 2 6> 2 = E Q a + 9.9 x lO^S" 2 , where 
S = R8 is the twisting-related arc length. Similar results 
are observed for 8 ss 60° as shown in Fig. [51 where we 
obtain E a = E™ - 8.2 x \Q- i R 2 {8 - 60) 2 . E™ is the 
value at 6 = 60°. 

In the above, we have presented the angle dependence 
for the interlayer potential of the twisting BLG. The re- 
mainder of this section will be devoted to explaining the 
origin for this angle dependence. Let's consider a single 
carbon atom on top of an infinite graphene sheet. The 
distance from this single atom to the graphene is fixed 
to be 3.478 A, which is the platform value in the inset 
of Fig. [3] We refer to such system as single atom (SA) 
model. The interlayer energy per atom in the SA model 
is calculated by Eg a = E tot /2, where E tot is the total in- 
terlayer energy between the single atom and the infinite 
graphene. Fig. [6] (a) shows the energy distribution of the 
SA model. The x and y axes are the xy position of the 
single atom. There arc translational and six-fold rota- 
tional symmetries in the energy distribution correspond- 
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FIG. 5: (Color online) The energy per atom versus twisting 
angle around 9 = 60°. For curves from top to bottom, the 
radius increases from 10 to 100 A. All curves can be well fitted 
to E a = E°+P(6- 60) 2 , where E° is the value at 6 = 60°. p 
is the fitting parameter. Inset shows the radius-dependence 
of p. P is fitted to P = -8.2 x 10~ 4 .R 2 . 



ing to the honeycomb structure of the graphene sheet. 
Fig. IH1 (b) shows a three-dimensional plot of the energy 
distribution within one hexagonal area. The hexagon is 
formed by six carbon atoms. When the single atom is 
on top of the hexagon corner (position A), E^a reaches 
the maximum value of E^a = Ea = —16.045 meV. 
When the single atom is on top of the hexagon cen- 
ter (position B), Eg a achieves the minimum value of 
E SA = E B = -29.810 meV. We note that the origin 
of the coordinate system here is set to the center of the 
hexagon in the bottom layer. 

For an arbitrary twisting angle, the positions of all 
atoms are arbitrarily distributed between position A and 
B. The energy per atom of such a system should be equiv- 
alent to the average of the energy in Fig. [6] (b) over the 
xy area within one hexagon: 



< E >- 



J E SA (x,y)dxdy 
J dxdy 



(2) 



We use a rectangular grid to partition the xy area, with 
N x and N y points in the x and y directions. The num- 
ber of total grid points is N x x N y . Fig. [5] (c) shows the 
convergence of the average with increasing grid points 
N x x N y . The average converges to a saturate value, 
which is exactly the energy per atom in Fig. [3] (the dashed 
red line). We propose a mapping between twisting BLG 
and the SA model: different twisting angle in the BLG 
corresponds to different grid type in the integration of 
Eq. 0). According to the Riemann-Lebesgue theorem^ 
the integration in Eq. ^ exists and does not depend on 
the grid type, because the integral function is bounded 
and smooth everywhere. From the mapping, this the- 
orem tells us that the energy per atom is the same for 
twisting BLG with twisting angle 9 in a large angle range. 
It explains why E a does not depend on the twisting angle 
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FIG. 6: (Color online) The registry-dependent potential be- 
tween a single atom and an infinite large graphene sheet. The 
single atom is on top of the graphene at a fixed distance 
z = 3.478 A. x and y axes in the figure are the other two 
coordinates of the single atom, (a) shows six-fold symme- 
try in the energy, due to the hexagon structure of graphene. 
(b) shows only one valley of the energy, (c) shows the conver- 
gence of the average of the energy over xy area with increasing 
grid points. The dashed line (red online) depicts the platform 
value of the energy per atom in Fig. [3] 



in a wide range. 

It should be note that the SA model can not be applied 
to explain the two boundary regions 9 « and 60° . It is 
because the SA model depends on the interlayer space, 
while the space is sensitive to the twisting angle in these 
two boundary regions around 9 = and 60°. As will be 
shown bellow, if the interlayer space is independent of 
the twisting angle, then the SA model will succeed for all 
twisting angles, including 9 w and 60°. 

Summary for this section. Using the registry- 
dependent potential, we have shown that the energy of 
the twisting BLG is insensitive to the twisting angle in 
a large angle range, which can be analyzed by the SA 
model. It illustrates that, from the energy point of view, 
there is no favorable twisting angles for the twisting BLG 
with 9 <S [10,50°]. Particularly, Moire pattern, i.e struc- 
ture with a commensurate twisting angle does not corre- 
spond to the energy minimum structure. 
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FIG. 7: (Color online) The interlayer potential energy calcu- 
lated from Lennard- Jones potential with different cut-off. The 
y-axis (-Ea) is the energy per atom. In this calculation, the 
radius of the top layer is 100 A. The twisting angle 8 = 21.8°. 
For a cut-off distance of r c = 100 A, the variation in the 
energy is on the order of 10 -5 meV. 
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FIG. 9: (Color online) The energy per atom versus twisting 
angle 9. The top inset shows the close-up of the middle region 
9 £ [20, 40]°, where curves for large radius of 500 and 1000 A 
are indistinguishable. Two bottom insets show the close-up 
of the boundary regions 9 « amd 60° . 



FIG. 8: (Color online) The saturation of the energy per atom 
with increasing radius for different twisting angles 9. (a) 
9 = 0, i.e AB-stacking BLG. (b) 9 = 60°, i.e AA-stacking 
BLG. The energy per atom in the AA-stacking BLG is radius 
independent, (c) 9 = 21.8°, i.e the first commensurable an- 
gle, (d) 9 = 23°. (e) 9 = 13.2°, i.e the second commensurable 
angle, (f) 9 = 11°. 



III. RESULTS AND DISCUSSION FOR 
LENNARD-JONES POTENTIAL 

Although Lcnnard- Jones potential is insufficient for 
the description of intcrlayer interaction in BLG, in this 
section, for the convenience of comparison, we present 
simulation results based on the Lennard- Jones poten- 
tial, V(r) = 4e((cr/r) 12 - (cr/r) 6 ), with e = 2.5 
meV and a = 3.37 A. The two parameters a and 
e are fitted to experimental values for the interlayer 
space and the phonon dispersion along FA direction in 
three-dimensional graphite4^ Usually, the cut-off in the 
Lennard- Jones potential is set to be around 10 A. Fig. [7] 
shows the energy per atom for a BLG with top layer 
radius as 100 A and the twisting angle 8 = 21.3°. It 
shows that the variation is still on the order of 0.1 meV 
for cut-off 10 A. This variation is too large for present 
study. For a cut-off 100 A, the variation in the E a is on 
the order of 10 -5 meV. In the following calculation, we 
use a cut-off 100 A for Lennard- Jones potential, so that 
our calculation precision is 10~ 5 meV. The required pre- 
cision is much higher for Lennard- Jones potential than 
the registry-dependent potential, because the Lennard- 
Jones potential gives much weaker twisting energy. With- 
out losing universality, we use the same interlayer space, 
3.35 A, for BLG with all twisting angles in this section. 



Fig. [8] shows the convergence of E a with increasing 
radius (i?) in BLG at various twisting angles. The twist- 
ing angles are 0.0, 60°, 21.8°, 23.0°, 13.2°, and 11.0° in 
panels from (a) to (f). Panels (d) and (f) correspond to 
two arbitrary twisting angles, while all other four pan- 
els correspond to commensurable angles. For all panels, 
the energy saturates at radius around 500 A. In panel 
(a), AB-stacking BLG, E a converges to a saturate value 
of -20.54613 meV. In panel (b), AA-stacking BLG, the 
saturate value is -20.11223 meV. The energy difference 
between these two stacking styles is 0.43390 meV, which 
is one order smaller than the value from the registry- 
dependent potential in the previous section. In all other 
four panels, E a converges to the same saturate values of 
-20.40857 meV. It becomes interesting that the energy 
in a BLG with commensurable angles shown in panels 
(c) and (e) are the same as that in BLG with arbitrary 
twisting angles shown in panels (d) and (f). It illustrates 
that the commensurable twisting angle does not corre- 
spond to the energy minimum state. Hence, from the 
energy point of view, the twisting angle observed in the 
experiment is not necessarily the commensurable angle. 
Instead, it can be an arbitrary angle. We have observed 
similar phenomenon in previous section, where the inter- 
layer interaction is described by the registry-dependent 
potential; so this finding does not depend on the poten- 
tial type. 

Fig. [S] shows that E a does not depend on the twisting 
angle in a wide angle range, except in the two bound- 
ary regions (9 « 0, 60°). For small radius, there is some 
obvious oscillation in the range of 8 £ [10, 50]°. With in- 
creasing radius, this oscillation amplitude decreases and 
becomes indistinguishable after R > 500 A. It indicates 
that the oscillation is actually due to the size effect, and is 
not related to the twisting angle. The top inset shows the 
close-up of the curve in the angle range of 9 £ [20,40]°, 
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FIG. 10: (Color online) The energy per atom versus twisting 
angle around 8 = 0. For curves Irom bottom to top, the radius 
increases as 20, 55, 150, 190, 250, 310, 400, 520, 665, 850, and 
1000 A. All curves can be well fitted to E a = E° + a8 2 , 
where is the value at 8 = 0. a is the fitting parameter. 
Inset shows the radius-dependence of a. a is fitted to a = 
3.9 x 1Q- 5 R 2 . 



FIG. 11: (Color online) The energy per atom versus twisting 
angle around 8 = 60°. For curves from top to bottom, the 
radius increases as 20, 55, 150, 190, 250, 310, 400, 520, 665, 
850, and 1000 A. All curves can be well fitted to E a = E° + 
13(8 - 60) 2 , where E° a is the value at 8 = 60°. /3 is the fitting 
parameter. Inset shows the radius-dependence of /3. /3 is fitted 
to p= -9.2 x 10~ 5 i? 2 . 



where the two results for large radius of 500 and 1000 A 
are indistinguishable. According to the energy minimum 
condition criteria, these results imply that the twisting 
angle can be an arbitrary value in this range. Two bot- 
tom insets show the close-up of the boundary regions 
9 « amd 60°. 

These two boundary regions are further zoomed in in 
Fig. rjTJ]and Fig. [TT] Fig. [TO] shows that E a can be fitted 
to E a = E® + a9 2 , where E° is the value at 9 = 0. a 
is the fitting parameter and is radius-dependent. Inset 
shows the radius-depcndcncc of a, which is fitted to a = 
3.9 x 10~ 5 i? 2 . Hence, around 9 = 0, we have E a = 
E° + 3.9 x 10- 5 R 2 9 2 = E° a + 3.9 x lO" 5 ^ 2 , where S = 
R9 is the twisting-related arc length. Similar results are 
observed for 9 sa 60° as shown in Fig.rjTJ where we obtain 
E a = El° - 9.2 x IQ- 5 R 2 (9 - 60) 2 . E™ is the value at 
9 = 60°. 

We are now applying the SA model to analyze these 
simulation results. In the SA model, the distance from 
the single atom to the graphene is fixed to be 3.35 A. 
We calculate the energy per atom in this SA model. 
Fig.[TOJ(a) shows the energy distribution of the SA model. 
Fig. [T2] (b) shows a three-dimensional plot of the en- 
ergy distribution within one hexagonal area. When the 
single atom is on top of the hexagon corner (position 
A), Esa reaches the maximum value of Eg a = Ea = 
—20.11223 meV. When the single atom is on top of the 
hexagon center (position B), £?sa achieves the minimum 
value of E SA = E B = -20.98000 meV. 

In the AB-stacking BLG, half of the atoms are located 
at position A and the other half atoms are located at 
position B. As a result, we get E a = (Ea + Eb) /2 = 
—20.54612 meV, which agrees quite well with the satu- 
rate value from Fig. [5] (a). In the AA-stacking BLG, all 



atoms are at position A, so E a = Ea = —20.11223 meV, 
which is exactly the same as the value from Fig. [5] (b). 
For an arbitrary twisting angle, the positions of all atoms 
are arbitrarily distributed between position A and B. The 
energy per atom of such a system should be equivalent 
to the average of the energy in Fig. [T^] (b) over the xy 
area within one hexagon following Eq. ©. Fig. Q~2] (c) 
shows the convergence of the average with increasing grid 
points N x x N y . The average converges to a saturate 
value, which is exactly the energy per atom in Fig. [3] 
(the dashed red line). 

We note that the interlayer space here has been kept 
the same for all twisting angles, so the SA model (with 
space 3.35 A) can be used for all twisting angles, in- 
cluding the two boundary regions 9 m and 60°. The 
energy curve in Fig. [12] (b) within the vicinity of the val- 
ley or the peak can be well described by functions E a = 
E min + 1.16r 2 - 0.44r 4 and E a = E max - 0.66r 2 + 0.53r 3 , 
respectively, where the distance r is with respect to the 
minimum point at the valley or the maximum point in 
the peak. 

For a twisting BLG, if the twisting angle 9 = 0, then 
half of the carbon atoms in the top layer are located at 
position B, i.e the minimum point of the valley; while the 
other half of atoms are at position A, i.e the maximum 
point of the peak. If the twisting angle is small but close 
to 0, i.e 9 w 0, half of the carbon atoms in the top layer 
deviate slightly from the valley while the other half of 
atoms deviate slightly from the maximum point of the 
peak. The deviation of each atom is s w r * 9, where r 
is the radial position of the atom. The edge atoms have 
the largest deviation S w R9, where R is the radius of 
the top graphene layer. 

It is clear that the energy per atom for the first half 
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suits (9.2 x 10~ 5 ) in Fig. [TT] Now it also becomes clear 
that the discrepancy in large radius limit between the fit- 
ting curve and the calculated data in Fig.fTOland Fig.fTTI 
is due to the nonlinear potential energy in higher order 
terms of (R9). 



FIG. 12: (Color online) The Lennard- Jones potential between 
a single atom and an infinite large graphene sheet. The single 
atom is on top of the graphene at a fixed distance z — c — 
3.35 A. (a) shows six-fold symmetry in the energy, due to the 
hexagon structure of graphene. (b) shows only one valley of 
the energy, (c) shows the convergence of the average of the 
energy over xy area with increasing grid points. The dashed 
line (red online) depicts the platform value of the energy per 
atom in Fig. [9] 

atoms can also be calculated by averaging of Eq. ([2]) over 
a circular area around the valley with radius S in the 
xy plane, yielding < E >= E min + 1.76 x 10~ 4 (R9) 2 - 
1.4 x 10~ 8 (i?6>) 4 , where 9 is in the unit of degree. For the 
other half atoms, the energy per atom can be obtained by 
doing average of Eq. ^ over a circular area around the 
peak with radius S in the xy plane. The integral result 
is < E >= El - 1.01 x \<T\my + 1.13 x 10- 6 (R6) 3 . 
As a result, the energy per atom for the twisting BLG 
with small twisting angle is < E >ss {E m [ n + E max )/2 + 
3.8 x 10~ 5 (R9) 2 , where the coefficient 3.8 x 10~ 5 agrees 
well with the results (3.9 x 1CT 5 ) in Fig. [TDJ Similarly, 
for twisting BLG with twisting angle 9 around 60°, the 
energy of all atoms deviate slightly from the maximum 
point of the peak, so the energy per atom from the in- 
tegral of Eq. © is < E >« £ max - 1.01 x \Q- A {R9f, 
with the coefficient —1.01 x 10~ 4 quite close to the re- 



IV. CONCLUSTION 

In conclusion, using both registry-dependent potential 
and Lennard- Jones potential, we have shown that the in- 
terlayer energy in a twisting BLG is independent of the 
twisting angle except in the boundary regions 9 ~ 0, 
60°. In these two boundaries, the interlayer energy is 
related to the square of the twisting-related arc length 
(S), E a = E® + cS 2 , with a constant coefficient c. These 
observations are well explained by averaging the energy 
distribution in the SA model, where a BLG with a par- 
ticular twisting angle is mapped to a special grid type 
in the partition of the integration area in the averaging. 
The energy distribution function in the SA model is so 
smooth that its integration docs not depend on the grid 
type, leading to the twisting angle independence of the 
interlayer potential in a twisting BLG. Our theoretical re- 
sults indicate that there is no preference for a commen- 
surable twisting angle in theBLG sample, according to 
the minimum energy condition criterion. These findings 
provide a possible explanation for the diverse twisting 
angles observed in the experiment. 
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